#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TSystem.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TEllipse.h"
#include "TClonesArray.h"

using namespace std;

void read(){

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  TFile *fread=new TFile("/scratch/shinichi/OutputFiles/75L_2Phase.root");
  TTree *tree=(TTree*)fread->Get("tree");
  
  TClonesArray *JPMTa=new TClonesArray("PMT");
  TClonesArray *MPMT1a=new TClonesArray("PMT");
  TClonesArray *MPMT2a=new TClonesArray("PMT");
  TClonesArray *MPMT3a=new TClonesArray("PMT");
  TClonesArray *MPMT4a=new TClonesArray("PMT");
  int nPhoton;

  tree->GetBranch("JPMTb")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT1b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT2b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT3b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT4b")->SetAutoDelete(kFALSE);
  tree->SetBranchAddress("JPMTb",&JPMTa);
  tree->SetBranchAddress("MPMT1b",&MPMT1a);
  tree->SetBranchAddress("MPMT2b",&MPMT2a);
  tree->SetBranchAddress("MPMT3b",&MPMT3a);
  tree->SetBranchAddress("MPMT4b",&MPMT4a);
  tree->SetBranchAddress("nPhoton",&nPhoton);

  double LowEdge = -500;
  double HighEdge = 100000;

  TH1F *hJPMTx=new TH1F("hJPMTx","hJPMTx",100,-150,150);
  TH1F *hJPMTy=new TH1F("hJPMTy","hJPMTy",100,-150,150);
  TH1F *hJPMTz=new TH1F("hJPMTz","hJPMTz",100,-650,-550);
  TH1F *hMPMTx=new TH1F("hMPMTx","hMPMTx",100,-150,150);
  TH1F *hMPMTy=new TH1F("hMPMTy","hMPMTy",100,-150,150);
  TH1F *hMPMTz=new TH1F("hMPMTz","hMPMTz",100,-650,-550);

  const int NEV=512;
  TH2F *hPMTxy[NEV];

  TH2F *hPosHit_All=new TH2F("hPosEff_All","hPosEff_All",30,-150,150,35,-650,-300);
  TH2F *hPosHit_J=new TH2F("hPosEff_J","hPosEff_J",300,-150,150,380,-680,-300);
  TH2F *hPosHit_M1=new TH2F("hPosEff_M1","hPosEff_M1",300,-150,150,380,-680,-300);
  TH2F *hPosHit_M2=new TH2F("hPosEff_M2","hPosEff_M2",300,-150,150,380,-680,-300);
  TH2F *hPosHit_M3=new TH2F("hPosEff_M3","hPosEff_M3",300,-150,150,380,-680,-300);
  TH2F *hPosHit_M4=new TH2F("hPosEff_M4","hPosEff_M4",300,-150,150,380,-680,-300);
  TH1F *hMeanAll=new TH1F("hMeanAll","hMeanAll",1000,0,1);
  TH1F *hMeanJ=new TH1F("hMeanJ","hMeanJ",1000,0,1);
  TH1F *hMeanM1=new TH1F("hMeanM1","hMeanM1",1000,0,1);
  TH1F *hMeanM2=new TH1F("hMeanM2","hMeanM2",1000,0,1);
  TH1F *hMeanM3=new TH1F("hMeanM3","hMeanM3",1000,0,1);
  TH1F *hMeanM4=new TH1F("hMeanM4","hMeanM4",1000,0,1);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;
  int cc=0;

  int mode = -1;
  if     ( nev==399 )   mode = 0;
  else if( nev==475 )   mode = 1;
  else if( nev==11286 ) mode = 2;
  else if( nev==17820 ) mode = 3;
  else{cout<<"!!! Wrong NEV = "<<nev<<" !!!"<<endl; return;}

  for(int iev=0;iev<nev;iev++){
    JPMTa->Clear();
    MPMT1a->Clear();
    MPMT2a->Clear();
    MPMT3a->Clear();
    MPMT4a->Clear();
    tree->GetEntry(iev);

    //hPMTxy[iev]=new TH2F(Form("hPMTxy_%d",iev),Form("hPMTxy_%d",iev),100,-150,150,100,-150,150);

    int njhit = JPMTa->GetEntriesFast();
    for(int ijhit=0;ijhit<njhit;ijhit++){
      PMT *jpmt = (PMT*)JPMTa->At(ijhit);
      hJPMTx->Fill( jpmt->xpos );
      hJPMTy->Fill( jpmt->ypos );
      hJPMTz->Fill( jpmt->zpos );
      //hPMTxy[iev]->Fill( jpmt->xpos,jpmt->ypos);
    }
    int nmhit1 = MPMT1a->GetEntriesFast();
    for(int imhit=0;imhit<nmhit1;imhit++){
      PMT *mpmt1 = (PMT*)MPMT1a->At(imhit);
      hMPMTx->Fill( mpmt1->xpos );
      hMPMTy->Fill( mpmt1->ypos );
      hMPMTz->Fill( mpmt1->zpos );
      //hPMTxy[iev]->Fill( mpmt1->xpos,mpmt1->ypos);
    }
    int nmhit2 = MPMT2a->GetEntriesFast();
    int nmhit3 = MPMT3a->GetEntriesFast();
    int nmhit4 = MPMT4a->GetEntriesFast();

    int iz,ir;
    int ix,iy;
    float xpos,ypos;
    float R = 0;
    if( mode==0 ){
      iz = iev / 21 % 19;
      ir = iev % 21;
    }
    if( mode==1 ){
      iz = iev / 25 % 19;
      ir = iev % 25;
    }
    if( mode==2 ){
      ix = iev % 22;
      iy = iev / 22 % 27;
      iz = iev / 22 / 27 % 19;
      xpos = -130.+ix*10.;
      ypos = -130.+iy*10.;
      R = sqrt( xpos*xpos+ypos*ypos );
    }
    if( mode==3 ){
      ix = iev % 22;
      iy = iev / 22 % 27;
      iz = iev / 22 / 27 % 30;
      xpos = -130.+ix*10.;
      ypos = -130.+iy*10.;
      R = sqrt( xpos*xpos+ypos*ypos );
    }

    //    float zpos = iz*10.-545.;
    float zpos;
    if( iz<21 )    zpos = iz*10.-565.;
    else           zpos = -359.+(iz-21)*1.;
    float rpos = ir*10.-120.;

    if( R>128 ) continue;
    hMeanAll->Fill( (njhit+nmhit1+nmhit2+nmhit3+nmhit4)/(float)nPhoton );
    hMeanJ  ->Fill( njhit/(float)nPhoton );
    hMeanM1 ->Fill( nmhit1/(float)nPhoton );
    hMeanM2 ->Fill( nmhit2/(float)nPhoton );
    hMeanM3 ->Fill( nmhit3/(float)nPhoton );
    hMeanM4 ->Fill( nmhit4/(float)nPhoton );
    cc++;

    if( ypos!=0 ) continue;
    if(xpos==0.&&zpos==-565) cout<<"N of Generated Photon = "<<nPhoton<<endl;

    hPosHit_All->Fill( xpos,zpos,(njhit+nmhit1+nmhit2+nmhit3+nmhit4)/(float)nPhoton );
    hPosHit_J  ->Fill( xpos,zpos, njhit/(float)nPhoton );
    hPosHit_M1 ->Fill( xpos,zpos,nmhit1/(float)nPhoton );
    hPosHit_M2 ->Fill( xpos,zpos,nmhit2/(float)nPhoton );
    hPosHit_M3 ->Fill( xpos,zpos,nmhit3/(float)nPhoton );
    hPosHit_M4 ->Fill( xpos,zpos,nmhit4/(float)nPhoton );

    
  }

  cout<<cc<<endl;
  cout<<" All Eff(Mean) = "<<100*hMeanAll->GetMean()<<" %"<<endl;
  cout<<"JPMT Eff(Mean) = "<<100*hMeanJ->GetMean()<<" %"<<endl;
  cout<<"MEG1 Eff(Mean) = "<<100*hMeanM1->GetMean()<<" %"<<endl;
  cout<<"MEG2 Eff(Mean) = "<<100*hMeanM2->GetMean()<<" %"<<endl;
  cout<<"MEG3 Eff(Mean) = "<<100*hMeanM3->GetMean()<<" %"<<endl;
  cout<<"MEG4 Eff(Mean) = "<<100*hMeanM4->GetMean()<<" %"<<endl;


  TCanvas *c1=new TCanvas("c1","c1",1000,800);

  TH1 *frame1=gPad->DrawFrame(0,0,0.03,3200);
  frame1->SetXTitle("Light Collection Efficiency");
  //  frame1->SetXTitle("N of Detected Photon (Er=300keV)");
  frame1->SetYTitle("NEV");

  hMeanJ->SetLineColor(1);
  hMeanJ->SetLineWidth(2.5);
  hMeanJ->Draw("same");
  hMeanM1->SetLineColor(2);
  hMeanM1->SetLineWidth(2.5);
  hMeanM1->Draw("same");
  hMeanM2->SetLineColor(4);
  hMeanM2->SetLineWidth(2.5);
  hMeanM2->Draw("same");

  TLegend *leg1=new TLegend(0.5,0.6,0.8,0.9);
  leg1->SetFillStyle(0);
  TLegendEntry *le=leg1->AddEntry( hMeanJ,"J-PMT","l");
  le->SetTextColor(1);
  le=leg1->AddEntry( hMeanM1,"Meg1 (Meg4)","l" );
  le->SetTextColor(2);
  le=leg1->AddEntry( hMeanM2,"Meg2 (Meg3)","l" );
  le->SetTextColor(4);
  leg1->Draw();

  /*
  TCanvas *c2=new TCanvas("c2","c2",1200,1000);
  c2->Divide(5,4);
  TEllipse *eJPMT=new TEllipse(0.,0.,32.,32.);
  TEllipse *eMPMT=new TEllipse(-2*40,0.,23.,23.);
  TCanvas *c3=new TCanvas("c3","c3",1200,1000);
  c3->Divide(5,4);

  for(int i=0;i<20;i++){
    c2->cd(i+1);
    hPMTxy[i]->Draw("colz");
    eJPMT->SetFillStyle(0);
    eMPMT->SetFillStyle(0);
    eJPMT->Draw();
    eMPMT->Draw();
    c3->cd(i+1);
    hJPMTt[i]->SetXTitle("time (ns)");
    hJPMTt[i]->SetYTitle("#P.E.");
    hJPMTt[i]->Draw();
  }
  */
  gStyle->SetOptStat(000000);

  double max = 0.015;
  TCanvas *c4=new TCanvas("c4","C4",800,600);
  //  c4->Divide(3,3);
  //  c4->cd(5);
  hPosHit_J->SetMaximum( max );
  hPosHit_J->SetXTitle("Y (mm)");
  hPosHit_J->SetYTitle("Z (mm)");
  hPosHit_J->Draw("colz");
  //  c4->cd(3);
  TCanvas *c5=new TCanvas("c5","C5",800,600);
  hPosHit_M1->SetMaximum( max );
  hPosHit_M1->SetXTitle("Y (mm)");
  hPosHit_M1->SetYTitle("Z (mm)");
  hPosHit_M1->Draw("colz");
  //  c4->cd(1);
  TCanvas *c6=new TCanvas("c6","C6",800,600);
  hPosHit_M2->SetMaximum( max );
  hPosHit_M2->SetXTitle("Y (mm)");
  hPosHit_M2->SetYTitle("Z (mm)");
  hPosHit_M2->Draw("colz");
  //  c4->cd(7);
  TCanvas *c7=new TCanvas("c7","C7",800,600);
  hPosHit_M3->SetMaximum( max );
  hPosHit_M3->SetXTitle("Y (mm)");
  hPosHit_M3->SetYTitle("Z (mm)");
  hPosHit_M3->Draw("colz");
  //  c4->cd(9);
  TCanvas *c8=new TCanvas("c8","C8",800,600);
  hPosHit_M4->SetMaximum( max );
  hPosHit_M4->SetXTitle("Y (mm)");
  hPosHit_M4->SetYTitle("Z (mm)");
  hPosHit_M4->Draw("colz");

}

